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Abstract. - The energy spectrum of the quantum Klein-Gordon (KG) lattice is computed 
numerically for model parameters relevant to optical phonon spectra. A pairing of phonon 
states is found when nonlinearity is significant, which agrees with other studies on different 
quantum lattice models [1, 2]. It results from the lattice anharmonicity, the magnitude of 
which is quantified by the binding energy of phonon bound states. Our work focuses on the 
case of weak anharmonicity, i.e., the phonon binding energy is weaker than the single phonon 
band width. We find that the phonon pairs dissociate at the center of the lattice Brillouin 
zone, whereas at the edge the binding energy remains comparable to the width of the single 
phonon band. Consequently, a weak nonlinearity is characterized by a pseudogap in the energy 
spectrum of two-phonon states. 



The phonon frequencies are some fundamental quantities that provide information about 
atomic interactions and local structure in crystals as well as in polymers or proteins. At the 
atomic scale, the phonon frequencies are formally derived from the quadratic expansion of the 
potential energy with respect to atomic displacements. Actually a more accurate description 
of the potential energy might be obtained by a higher order expansion which involves some 
non-quadratic terms. Those terms are usually referred to as nonlinear terms because they 
yield some forces that are not proportional to displacements. In the early seventies, the non- 
quadratic contribution to the energy has been quantified in molecular crystals such as CO2, 
N2O and OCS (see Ref. [3] and Refs. therein). The nonlinearity was then identified by 
some anharmonic peaks in the infrared spectrum. In solid i^, similar infrared resonances 
were earlier interpreted [4] as a signature of phonons pairing. Currently the list of materials 
in which phonon bound states occur is still growing [5,6]. In addition, some controversial 
cases, as for instance acetanilide [7] or the spectral branch A^LO) of diamond [8,9] can be 
mentioned as still open questions. In the present study, we point out a possible difficulty for 
the vibration spectrum interpretation which could occur when nonlinearity is not sufficiently 
strong to separate the anharmonic resonances from the harmonic energy regions. Then, the 
hybridization between bound and unbound phonons is shown to imply a pseudogap in the 
lattice energy spectrum. 

The anharmonicity of molecular crystals has been theoretically studied by V.M. Agra- 
novich [1] , who proposed a qualitative approach where some boson quasi-particles model the 
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atom vibrations. A Hubbard onsite interaction between boson pairs was introduced to sim- 
ulate the nonlinear behavior of the lattice (for recent studies see Ref. [2, 10]). This effective 
boson model displayed the concepts of phonon bound states and more specifically of biphonon. 
The drawback of the model is that the energy terms that do not conserve the boson number 
are neglected, although those terms stem from the potential energy of atoms (or molecules). 
A rigorous approach to models which do not have a conserved boson number was proposed in 
Ref. [11]. The present paper can be considered as numerical support for that theory. Here, a 
numerical method is introduced for computing the KG energy spectrum with different forms 
of nonlinearity, i.e., different orders in the potential energy expansion. Our technique is based 
on the construction of a nonlinear phonon basis which allows to treat lattices whatever the 
magnitude of nonlinearity is, in contrast to the linear phonon basis [12], derived from the 
harmonic approximation. The attractive advantage of our method is that the maximum size 
of the lattices that can be studied is large enough to approach the infinite system features, i.e., 
the energy spectrum shows no boundary effects and it is nearly continuous in the reciprocal 
space. 

The KG model may be introduced as follows. At the node i of a translational invariant 
d-dimensional lattice, the internal mode Xi of a unit cell (a molecule, for instance) evolves in 
the local potential V with an effective mass to, each of the xi being coupled to the nearest 
neighbors' modes by a quadratic coupling, parametrized by c. Then, the Hamiltonian of that 
system reads: 

H = T,^ + v ^- c E ( x *- x i) 2 ]- (!) 

i j=<i> 

The subscript j =< i > denotes the first neighbors of site i. The operator Xi and its conjugate 
momentum pi satisfy [xi,pi] — ih. For small amplitudes of Xi, the Taylor expansion of V gives 
V(xi) = a,2xf + a^xf + a^xf . It has been limited to the fourth order but higher order terms can 
be taken into account with no difficulties in calculations that follow. Higher order coupling 
terms have been found not to modify qualitatively our results. Introducing the dimensionlcss 
coefficients and operators: 

A a I ^ A - a ^ C - 4 ° 

Pi = Pi/VmhH, and Xi = Xi^mfl/h (2) 

where the frequency = \/2(a 2 — {2d)c)/m characterizes the d-dimensional oscillator net- 
works, the Hamiltonian reads: 

H = %nY,[^ + ^+A i Xf + A l Xt + ^X i E Xj]- (3) 

i j = <i> 

Note that the coupling coefficient c contributes to Q. The width of an optical phonon band is 
physically a few percent of the elementary excitation energy (w Ml). We restrict our study to 
A 4 > to ensure that the local potential V is positive when the onsite displacement increases 
arbitrarily. In contrast, A 3 can take positive or negative values provided V is a monotonic 
function. 

Our computation is now described. The first step is to calculate eigenstates of the onsite 
Hamiltonian hi — Pf /2 + Xf /2 + A$Xf + A 4 Xf, which is performed by projecting hi over the 
usual Einstein basis. The diagonalization of the corresponding matrix is realized numerically 
with great accuracy. The eigenstates of hi are arranged by increasing order of eigenvalues, 
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so the a th eigenstate is denoted <j> at i and its eigenvalue is 7(a). Considering the Hamiltonian 
Hq = HQ ^ i hi, the Hq eigenstates can be written as products of onsite states Tli4>ai.i- Taking 
advantage of the lattice translational degeneracy, some Bloch waves are constructed from those 
products as follows: 



where ao is the lattice parameter, q is the wave vector, the subscript [II;a;] identifies different 
types of products and A[n iai ] ensures the normalization. Expanding the Hamiltonian H in 
the Bloch waves is done analytically and gives a matrix B(q). Since H does not hybridize 
Bloch waves with different q, the diagonalization of B(q) can be performed independently for 
each q. To that aim, we used exact numerical methods (LAPACK or Numerical Recipes) on 
a simple desktop PC. Then, the Schrodinger equation is solved with an error which shrinks 
to zero exponentially by increasing the cutoff over the H eigenstates. The cutoff is fixed by 
5^ &i < N cut . The number of Bloch waves involved in diagonalization is 3052 for a lattice 
size S = 33. The accuracy of our calculations has been tested both on the anharmonic atomic 
chain 5 = 4 [12] and the harmonic chain for which the eigenvalue problem is solved analytically 
by a spatial Fourier transform of Eq. ^ For these two comparisons, very good agreements 
have been found. For the latter case, the results are reported in Fig. ^ Some distinct cutoff 
values N cu t have been tested and clearly our eigenvalue evaluation converges to the exact 
values as the cutoff increases. For two-phonon excitations, the error is less than 0.1% when 
N C ut — 4. Then the accuracy is even better for the phonon band which fits the dispersion law: 
fiVL^/l + 2C cos(qa n ). 

Our results are first presented for a ID lattice. In Figs. |2Ja-c), for a finite size S = 33, 
the eigen-energies (circle symbols) are plotted versus the wave momentum q. They contribute 
to the distinct branches of the spectrum. When the non-quadratic part of the energy is 
negligible, the eigen-spectrum is composed of the fundamental optical branch and the linear 
superpositions of phonon states. Increasing gradually A4, fixing C = 0.05 (relevant for optical 
phonons, see Fig. Ufa)) and A3 = 0, makes a branch split from the top of the two-phonon 
band. Then, the onsite potential V is said to harden. Subsequently, if A4 is fixed to a positive 
arbitrary value and A3 is increased in absolute value, the isolated branch backs into the two- 
phonon band and eventually splits from the bottom for large enough \A$\. The onsite potential 
V is now said to soften. By analogy with the biphonon theory [1] , the splitting branch (labelled 
by {2} in Figs. I2b-c)) is identified as being the energy of biphonon states. Measuring the 
energy of the biphonon with reference to the unbound two-phonon states for same momentum 
q, a biphonon binding energy can be defined. The absolute value minimum, with respect to q 
of the biphonon binding energy is called the biphonon energy gap. When the non-quadratic 
contribution is not sufficient to open the biphonon gap, the gap is found to close at the center 
of the lattice Brillouin zone (BZ) (Fig|5fc)). On the other hand, at the edge of BZ the binding 
energy may remain comparable to the width of the single phonon branch. Thus, a pseudogap 
occurs when the non-quadratic energy has the same magnitude as inter-site coupling. Similar 
pseudogaps have been found in other quantum lattices [13,14]. In the pseudogap regime, 
biphonon excitations exist only at the edge of BZ while they dissociate into unbound phonon 
pairs at center of the zone. In Figs. |5Jb-c), it is noteworthy that the hybridization between 
the biphonon states and the two-phonon states modifies substantially the biphonon branch 
curvature and consequently the quantum mobility of the biphonon. Such a feature can not 
be analyzed accurately within the effective Hubbard model for bosons since it neglects some 
energy terms that do not conserve the boson number. The same remark holds for the single 
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eigenvalue rank 

Fig. 1 - For a dimensionless coupling C = 0.05, comparison between the exact calculation (thick solid 
line) and our evaluations (detailed in the text) of the eigenvalues in a ID harmonic lattice composed 
of S = 33 oscillators. Our numerics are performed for 3 distinct cutoffs N cu t = 2 (dashed line, 
triangles), N cut = 3 (dot-dashed line, squares) and N cu t = 4 (thin solid line, circles). The spectrum is 
measured with respect to the groundstate and is arranged in increasing order. The eigenstate number 
is reported on the X-axis. The Y-axis unit is HQ. 

phonon dispersion law which takes the form (1 + Ccos(qao)) (with arbitrary energy units) 
in the effective boson model instead of the well-known y/l + 2Ccos(qao), calculated in the 
harmonic approximation for similar parameters. The two formulas diverge with increasing 
C, i.e. when the quantum hybridization is amplified. This disagreement is due to the boson 
product operators such as (af +1 af) that are neglected in the effective boson model, in contrast 
to the KG model and its harmonic approximation. 

Other computations have been performed with different model parameters to complete our 
results. The curvature of the biphonon branch depends on the sign of C provided the biphonon 
gap is large. Otherwise, the hybridization between bound and unbound phonons tends to 
impose the same curvature as the two-phonon band(Fig|2Ic)). The biphonon band curvature 
thus results from the interplay between the biphonon tunneling and the hybridization with 
unbound phonons. In Fig^b), the triphonon energy branch (labelled by {3}) which splits 
from the three-phonon band (not plotted) can also be noted with a similar behavior as the 
biphonon branch. For the same nonlinear coefficients as in FigsEJa-b), varying artificially C 
to the trivial case C = (see Fig.^J demonstrates that the isolated branches, {1}, {2} and {3} 
coincide with the energies of Bloch waves -B[ Q .o,...i with a single onsite excitation a — 1,2 and 
3, respectively. The branch tag corresponds to the excitation level {a} of the coinciding Bloch 
waves at C — 0. The reason for the isolated branches is the anharmonicity of hi eigenvalues 
j(a > 1) which implies some energy gaps at C = that persist when C ^ 0. In Fig. [21 we 
note the bands associated to unbound phonon pairs, to unbound three-phonon states and to 
unbound states composed of a biphonon and a single phonon. These states are labelled by 
multiple tags {a, (3, ...} that correspond to Bloch waves -B[a,...„8...] at C = 0. 
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Fig. 2 - Energy spectrum of a chain composed of S = 33 unit cells: (a) the optical phonon branch 
for model parameters C = 0.05, A3 = 0.13 and A4 = 0.01; (b) the two-phonon energy region for 
same parameters as (a); (c) same as (b) but A3 — 0.1. The biphonon branch is labelled by {2}, the 
unbound two-phonon bands by {11} and the triphonon branch by {3}. The Y-axis unit is TiQ. and 
its zero is groundstate energy. The wave vector q is reported on the X-axis, whose unit is (7r/ao). It 
ranges from the center to the edge of the first lattice Brillouin zone. 



Since our calculations are performed with large enough lattices to avoid boundary effects, 
there is a perfect agreement between the energy spectrum computed with S — 33 (Fig.[2a-b)) 
and with S = 19 (Fig.|3J for the same parameters. In Fig.|3 the band widths of phonon bound 
states increase with C much slower than for unbound phonons. The branches of biphonon 
{a — 2} and triphonon {a = 3} are found to merge with the unbound phonon bands above a 
certain threshold C a . For C < C a , the a th branch is separated from the rest of the spectrum 
by gaps that open over the whole BZ whereas around C ~ C a , only a pseudogap separates 
partially that branch from a certain unbound phonon band. The C a threshold depends on 
both A3 and A4 and it is different for each phonon bound state. When C is much larger 
than Ci , the biphonon gap is closed then only gaps of high order phonon bound states (such 
as triphonon) may be opened in the excitation spectrum. The results reported in Fig. [3] are 
qualitatively similar to the Raman analysis [15] of the high pressure molecular solid H2 which 
shows a pressure-induced bound-unbound transition of the so called bivibron, around 25 GPa 
(compare Fig. 10 in Ref. [15] and the two-phonon energy region in Fig.[3J). In the framework of 
our model, the pressure variation of experiments [15] can be simulated by a change of coupling 
integral C due to the fact that neighboring molecules are moved closer together because of 
the external pressure. The increase of C induces a bound-unbound transition of the biphonon 
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Fig. 3 - Plot of the energy spectrum versus dimensionless coupling C, for a ID quantum chain 
composed of S — 19 cells with model parameters A$ = 0.13 and A4 = 0.01. The energy unit is same 
as in Fig. [5] Labels are explained in the text. 



at C = Ci- The pressure variation of C as well as of other KG model parameters could 
inform us on how atomic potentials depend on the inter-atomic distances (that are known 
from diffraction measures). A biphonon is also exhibited by the internal stretching [5] of CO 
molecules adsorbed on the surface i?u(lll). There, the coverage of the surface can also be 
thought to change the coupling integral C, and a bound-unbound transition could thus be 
expected at high surface coverage. Yet, that assumption is seemingly inapplicable according to 
a recent theoretical work [10] which invoked a coverage defendant anharmonicity which would 
originate in a chemical modification of the intramolecular CO potential due to surrounding 
molecules. 

Finally, for a 2D lattice, a biphonon pseudogap is found to open around q — 1/2 [11] (see 
Pig. 01 . As we chose a set of nonlinear parameters such that the local potential V hardens, the 
biphonon energy is higher than for two-phonon instead of lower when V softens (see Fig. 
Aside from that point, the binding energy vanishes at center of BZ whereas the width of 
pseudogap at the edge of BZ has same order as the phonon branch width. That likeness 
between ID and 2D spectra shows that the pseudogap is a consistent property of lattices 
where nonlinearity is comparable to the inter-site coupling, whatever the dimensions are. For 
a 3D lattice, the biphonon pseudogap is expected to open around q = 1/2[111]. 

To summarize, the energy spectrum of the KG lattice has been considered with numerics. 
The KG Hamiltonian provides a realistic modelling of both quantum hybridization and pairing 
of phonon states. Through our computing method, the spectral features of phonon bound 
states can be related unambiguously to non-quadratic terms of the potential energy. We 
predict that when those terms are not strong enough to separate the spectral resonances 
of bound and unbound phonon states, the pseudogaps of the biphonon, triphonon or even 
higher order of phonon bound states, are potential signatures of nonlinearity. The biphonon 
pseudogap has been proved to open systematically at the edge of the first Brillouin zone. 
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Fig. 4 - Two-phonon spectrum for a 2D square lattice composed of S — 13 x 13 cells and for the 
model parameters C = 0.025, ^3 = and A4 = 0.025. Energy unit and tags are same as in Fig. [3] 
The color picture (easier to read) is available online. 
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